Symbolic Computation for Analysing Problems

TL;DR, Ditch pen/paper derivations and Excel monstrosities and use symbolic computation instead. Download MATLAB and this file, and have a play around.

Table of Contents

1. Introduction
2. Tutorial
   • Setup
   • Manipulation
   • Substitution
   • Calculation
   • Comparing Calculations
   • Interactivity/Flexibility
   • Optimisation
3. Documentation
4. Pitfalls
5. Tool Comparison
6. Other Symbolic Libraries
7. Conclusions

Introduction

Prior to discovering symbolic computation was a thing a few months ago, when making models, I would manipulate the relevant equation from a formula book until I get the form I need. I’d then make the decision on whether to use a calculator, Excel or MATLAB, depending on the complexity of the equation. Usually I would pick Excel.

As the scope of the analysis creeps, the spreadsheet usually becomes too difficult to use. It becomes even harder to decipher if you’re not the one who created the sheet…

Clearly there’s room for improvement…

Tutorial

For the tutorial, we’ll:

• Start with a constant acceleration equation of motion
• Manipulate it until it’s in the correct form
• Determine how long a ball will stay in the air if thrown directly upwards
• Convert the analysis from a 1D problem to a 2D problem
• See how the Initial throw angle affects the trajectory interactively
• Optimise the angle at which a ball should be thrown to cover the furthest distance

All that in ~60 lines of self-explanatory code touching neither pen nor paper. I’ve opted to do this in MATLAB for reasons explained later, using a trial version of the Symbolic Math toolbox.

If you don’t have time to follow the tutorial, you can download the file and view in MATLAB. You’ll still be able to view the Live Script exactly as shown in the GIF below, even if you don’t have the toolbox installed. You won’t be able to run the code if you make any edits unfortunately until you add the toolbox. Adding a trial version of the toolbox is as easy as clicking on the Add Ons button, and searching for Symbolic Math.

Setup

Once you’ve created a new Live Script, we create our algebraic variables with the command syms. These are called symbolic variables, and don’t have a numeric value. By default, MATLAB assumes all variables are complex unless you include assumptions such as “real” or “positive”.

>> syms s a t real

We can also initialise the symbols as a matrix of symbols (e.g. to signify those variables in different states or as a vector)

>> syms v [1 2] real

To create an equation, we write the following:

>> eqVUAT = v(2) == v(1) +a*t

That line of code states the second element of the matrix v (i.e. v(2)) is equal (==) to the first element of v (i.e. v(1) ) plus the product of acceleration and time i.e. (a*t). We then store that relationship in the variable eqVUAT using =.

Hit “Run” or “Ctrl + Enter” to output a beautifully typeset equation. Subscripts, superscripts, accents and Greek symbols are all supported.

eqVUAT =

If this would have been your final equation, and you’re ready to start calculating, skip to here.

Manipulation

Currently, that equation isn’t particularly useful for us, as we don’t know what v₂ or t is. We can integrate it to get a different form.

>> eqSUAT = s == int(rhs(eqVUAT),t,0,t)

eqSUAT =

In the above line of code we:

1. Use rhs() to extract the right hand side of the eqVUAT equation
2. Use int() to integrate the right side with respect to t, between the limits t = 0 and t = t
3. Use == to set the integral equal to the displacement s
4. Use = to store that equation in eqSUAT

Now we can rearrange to get t using isolate().

>> eqSUAT = isolate(eqSUAT,t)

eqSUAT =

N.B. Be aware that if there are multiple solutions (e.g. sqrt(x) and arcsin(x) ), it will take the simplest solution. Use solve() to return more solutions and pick the one you want. I’d recommend using it in a function as follows so that solve() returns an equation rather than just the right hand side.

function eqOut = isolate2(eqIn,var,solutionnumber)
    sol = solve(eqIn,var);
    eqOut = var == sol(solutionnumber);
end

Substitution

Maybe we’ve suddenly changed our minds, and now want to define s as the difference in start and end displacement, i.e. s = s₂ – s₁. Not a problem, we can do that using subs(). This function finds all instances of a specified variable/expression and replaces it with the variable/expression what you want.

>> syms s_1 s_2 real
>> eqSUAT = subs(eqSUAT,s,s_2-s_1)

Alternatively, if you’ve created that relationship already in another equation, you can substitute using lhs()and rhs() to extract the left and right hand sides of the equation.

>> syms s_1 s_2 real
>> eqS1S2 = s == s_2 - s_1;
>> eqSUAT = subs(eqSUAT,lhs(eqS1S2),rhs(eqS1S2))

eqSUAT =

Calculation

Now we can substitute numerical values. To do that, we create an object vars and store the value for each variable. Each value can be multiplied by u.unit to add pretty much any unit.

>> u = symunit;

>> vars.v1 = 30*u.mph;
>> vars.a = -9.81*u.m/u.s^2;
>> vars.s_1 = 0*u.m;
>> vars.s_2 = 0*u.m;

We then substitute in the values using the same subs() command as before.

>> subs(eqSUAT, vars)

ans =

Obviously, the above isn’t very useful. We need to follow a few steps to make it useful.

>> vpa(simplify(unitConvert(subs(eqSUAT,vars),"SI")),3)

ans =

1. Use subs() to substitute the values into the equation
2. Use unitConvert() to convert all the units to a specific unit system (“SI”) or unit (“u.s”)
3. Use simplify()simplify terms like 100*(30+sqrt(900)) to 6000
4. Use vpa() to switch to floating point arithmetic (e.g sqrt(2) to 1.414; 981/100 to 9.810)

Note that if you substitute a numeric value for t, MATLAB will check whether the left and right hand side match and output true or false. It will also output false if the solution is inconsistent with the assumptions, for example:

>> syms x y
>> assumeAlso(y,'real')
>> eqXY = y == sqrt(x);
>> simplify(subs(eqXY,x,-1))

ans=
FALSE

Comparing Multiple Calculations

Want to compare results from multiple different values quickly? No need to create an entire new array, just reassign the variable to include a list of values instead of a single value.

>> vars.v1 = [30,20]*u.mph;
>> vars.s_1 = [0*u.m,2*u.ft];

>> vpa(simplify(unitConvert(subs(eqSUAT,vars),"SI")),3)

ans =

Interactivity/Flexibility

I’ve now decided that I want to make the analysis 2D. I’ve done this by:

1. Using syms to initialise the new horizontal and vertical components.
2. Copying over the rearranged 1D equation to a new variable eqSUAT2D, in case you want to keep the old 1D equation for later
3. Using subs()to replace each variable (except for time) with a two variables corresponding to the horizontal and vertical components. e.g.:
   1. {s_2} is substituted by {s_2_h, s_2_v}
   2. {a} is substituted by {a_h, a_v}
   3. {v(1)} is substituted by {v(1)*cos(theta), v(1)*sin(theta)}

Once substituted, the output will be a 1×2 array of equations, corresponding to the horizontal and vertical equations of motion.

>> syms s_1_h s_1_v s_2_h s_2_v a_h a_v real
>> syms theta positive real
>> assumeAlso(a_v<0)
>> assumeAlso(a_h<=0)
>> assumeAlso(theta<pi/2)

>> eqSUAT2D = isolate(eqSUAT,s_2);
>> eqSUAT2D = subs(eqSUAT2D,{s_2, s_1,v(1),a},{[s_2_h, s_2_v],[s_1_h,s_1_v], 
   v(1)*cos(theta),v(1)*sin(theta)],[a_h,a_v]})

eqSUAT2D =

The above equations are in a parametric form that we can directly plot using fplot() between t = 0 and t = 5. Unfortunately, the fplot() function doesn’t accept units, so we leave them out entirely from the start or strip them out using separateUnits().

>> vars.a_h = 0*u.m/(u.s)^2;
>> vars.a_v = -9.81*u.m/(u.s)^2;
>> vars.v1 = 13.41*u.m/u.s;
>> vars.s_1_h = 0*u.m;
>> vars.s_1_v = 2*u.m;
>> vars.s_2_v = 0*u.m;
>> vars.s_2_h = s_2_h*u.m;
>> vars.theta = 45*u.deg;
>> vars.t = t*u.s;
>> vars.s1 = s1;

>> eqSUAT2Dsubs = separateUnits(unitConvert(subs(eqSUAT2D,vars),"SI"));

>> fplot(rhs(eqSUAT2Dsubs(1)),rhs(eqSUAT2Dsubs(2)),[0,5])
>> xlabel('Horizontal Distance (m)')
>> ylabel('Vertical Distance (m)')
>> xlim([0,20])
>> ylim([0,10])

To make this plot more interactive, highlight a number you’d like to modify and insert a slider:

Optimisation

If we wanted to figure out the optimum angle at which the ball should be thrown to cover the most distance, we determine when ds2h/dθ = 0.

>> eqSUAT2D

eqSUAT2D =

However, the time t when the ball lands is a function of θ, and hence needs to be substituted

>> eqTime = isolate2(eqSUAT2D(2),t,2)

eqTime =

>> eqDistanceTravelled = subs(eqSUAT2D(1),lhs(eqTime),rhs(eqTime))

eqDistanceTravelled =

Now we have the horizontal displacement in terms of constants and θ, we can differentiate the right hand side of the equation using diff() to get ds/dθ.

>> dfDistanceTravelled(theta) = simplify(diff(rhs(eqDistanceTravelled),theta))

dfDistanceTravelled(theta) =

MATLAB is unable to solve ds/dθ = 0 due to the complexity of the above equation. To solve it, we have a couple of options.

Numerically
We can solve it numerically, by substituting the actual values of all the variables except for θ, and use vpasolve(), giving a solution range of between 0 and pi/2.

>> vars.theta = theta; %ensuring theta is substituted with theta
>> dfDistanceTravelledsubs(theta) = subs(dfDistanceTravelled, vars)
>> vpasolve(dfDistanceTravelledsubs==0,theta,[0,pi/2])

ans = 
0.7361327600173396501181

which corresponds to 42.17°.

Analytically
Alternatively, we can solve it for a special case where the horizontal acceleration, start and end vertical displacements are zero (a_h, s_1_v, s_2_v = 0)

>> spec.a_h = 0;
>> spec.s_1_v = 0;
>> spec.s_2_v = 0;
>> dfDistanceTravelledSpecial(theta) = simplify(subs(dfDistanceTravelled,spec))

dfDistanceTravelledSpecial(theta) =

>> isolate(dfDistanceTravelledSpecial==0,theta)

ans =

Documentation

With Notebook style environments, we aren’t limited to documenting code using comments. We can write fully formatted text, which allows us to provide a much clearer narrative. The options include:

• Headings, sub headings, sub sub headings etc.
• Tables of Contents
• LaTeX equations in line with text
• Hyperlinks (internal and external)
• Images (Paste directly or through file dialog)
• Hiding all code, leaving just formatted text, and the outputs (e.g. equations and graphs)

In addition, we can split code up into different sections, so you can rerun sections with different parameters without running the entire script to speed things up a bit.

Pitfalls

Verification: The ease of use could lead you into a false sense of security…YOU WILL STILL HAVE TO VERIFY WHAT YOU DERIVE. However, this should be easier to verify as the equations are much easier to read, and interact with.

Cost: I used the 30 day trial of the symbolic math toolbox as CC doesn’t have a license for this toolbox yet. Prior to trying out the MATLAB version, I used Sympy, Jupyter Labs and Pint (a unit systems library) to the same effect though it’s not as cohesive as the MATLAB environment is, and a bit harder to share and setup. It’s all free though.

Setup Time: The setup time is quite fast, but possibly not as fast as Excel or Smath, but that doesn’t account for how Excel and Smath require you to solve things by pen/paper first.

Runtime: Solving things symbolically does take slightly more time. For instance, the running the entire script above takes ~3 seconds, and an additional second or so to render plots and equations whereas the equivalent floating point arithmetic script would take <0.5s. By splitting up the code into sections, this is significantly minimised however.

Tool Comparison

Here’s a subjective comparison of the different tools. Clearly the different tools do have their own strengths and you’ll have to decide yourself which one suits your use case best. For example, I’m uncertain how a tolerance analysis spreadsheet could be done in a notebook/Live Script environment well.

**You can still use non-symbolic, numeric arithmetic code in your Live Scripts/Notebooks, for all the speed benefits, just make sure to section correctly so you don’t run everything.

Initial Setup Speed: How easy is it to go from a formula to an answer?
Interactivity: How easy is it to modify parameters for a quick play around?
Flexibility: A kinda handwavy term for how easy it is for you to rework your analysis if you decide something else is important.
Optimisation: How easy is it to find the optimum parameters for your design?
Computation Time: How long does it compute an answer given an equation?
Readability: How easy is it for someone else to look at your analysis and figure out what’s going on? How easy is it for you to debug if something is wrong?

It is quite expensive though…

Conclusions

Since discovering this was a thing, I’ve barely used Excel as an analysis tool (except for the simplest use cases where a calculator would also do the job). In addition, because it’s so easy and quick to do, I’ve started analyses I would otherwise not have done at all due to time constraints…
There’s still a fair amount that the toolbox can do, that I didn’t cover:

• Matrix manipulation
• Expanding and factorising expressions, partial fraction decomposition,
• Manipulating trigonometric/hyperbolic/exponential identities
• Dimensional analysis
• Differential equations
• Solving systems of linear equations
• Symbolic Functions (f(x) etc.)

Here are some links to other scripts I’ve written:

• Power Transmission through an Interference Fit
• Analysing forces in a over-centre linkage
• Determining the critical speed to cause a ventilator to topple
• Optimising surface area of a box